############ A Second Experiment: Randomization in Only the Consent Form (Section J)

#### Libraries
library(sandwich)
library(lmtest)
library(ggplot2)
library(digest)
library(xtable)

# colour scheme
theme_bw1 <- function(base_size = 16, base_family = "") {  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text      
      axis.ticks =        element_line(colour = "grey75"),
      axis.title.x = element_text(size=base_size, face="bold"),
      axis.title.y =      element_text(size = base_size,angle=90, face="bold"),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}

#### Load working data
researcher.data.full <- read.csv("data_investigator_characteristics_second_experiment.csv")

### Was randomization successful (failed if name == Ebony Robinson)
researcher.data.full$random.success <- ifelse(researcher.data.full$researcherName == "Ebony Robinson",0, 1) ## No zeroes, so that's good.

### Convert names to treatment labels
blackwomen <- c('Tanisha Dorsey', 'Deja Gaines', 'Deja Banks', 'Latoya Dorsey', 'Aisha Rivers', 'Keisha Williams','Tanisha Rivers',
                'Keisha Joseph', 'Shanice Williams', 'Aisha Mosley')
whitemen <- c('Cody Novak','Neil Kelly','Connor Sullivan','Greg Sullivan','Brett Schwartz','Brad Walsh','Connor McCarthy',
              'Matthew Schmitt','Greg Murphy','Todd McCarthy')
blackmen <- c('Tyrone Booker','Rasheed Washington','Reginald Jones','Darnell Washington','Darnell Banks','Kareem Dorsey','Rasheed Jackson',
              'Darnell Booker','Rasheed Dorsey','Tyrone Gaines')
whitewomen <- c('Molly Koch','Molly Baker','Carrie Haas','Molly Sullivan','Abigail Koch','Kristen Murphy','Carrie McCarthy',
                'Molly Murphy','Abigail Ryan','Carrie Ryan')

# Treatments
researcher.data.full$treatment.whitewoman <- ifelse(researcher.data.full$researcherName %in% whitewomen, 1, 0)
researcher.data.full$treatment.blackwoman <- ifelse(researcher.data.full$researcherName %in% blackwomen, 1, 0)
researcher.data.full$treatment.whiteman <- ifelse(researcher.data.full$researcherName %in% whitemen, 1, 0)
researcher.data.full$treatment.blackman <- ifelse(researcher.data.full$researcherName %in% blackmen, 1, 0)

# Woman/Black treatments
researcher.data.full$treatment.woman <- ifelse(researcher.data.full$treatment.whitewoman == 1 | researcher.data.full$treatment.blackwoman == 1, 1, 0)
researcher.data.full$treatment.black <- ifelse(researcher.data.full$treatment.blackwoman == 1 | researcher.data.full$treatment.blackman == 1, 1, 0)

#### Did the respondent finish the survey?
researcher.data.full$NoFinish <- abs(1 -  researcher.data.full$Finished)


#############################
### Treatment effects

################################
#### Subset data to units that finished
researcher.data <- subset(researcher.data.full, Finished == 1)

#### All -99 are NAs
researcher.data[researcher.data==-99] <- NA

###################################################
### Pass attention checks?

researcher.data$ac1.pass <- ifelse(researcher.data$ac1.huffpo == 1 & researcher.data$ac1.Reuters == 1, 1, 0)
researcher.data$ac2.pass <- ifelse(researcher.data$ac2.Energy == 1 & researcher.data$ac2.GlobalTrade == 1, 1, 0)

researcher.data$ac1.pass[is.na(researcher.data$ac1.pass)] <- 0
researcher.data$ac2.pass[is.na(researcher.data$ac2.pass)] <- 0

researcher.data$ac.pass.both <- ifelse(researcher.data$ac1.pass== 1 & researcher.data$ac2.pass == 1, 1, 0)

####################################################
### Background covariates

researcher.data$gender.male <- ifelse(researcher.data$gender == 1, 1, 0)
researcher.data$gender.female <- ifelse(researcher.data$gender == 2, 1, 0)

researcher.data$democrat <- ifelse(researcher.data$party == 1, 1, 0)
researcher.data$republican <- ifelse(researcher.data$party ==2, 1, 0)
researcher.data$independent <- ifelse(researcher.data$party == 3, 1, 0)

researcher.data$voteyes <- ifelse(researcher.data$turnout == 3, 1, 0)
researcher.data$voteno <- ifelse(researcher.data$turnout == 2, 1, 0)
researcher.data$votedk <- ifelse(researcher.data$turnout == 1, 1, 0)

################################################3

###### Racial_resentment Scale
# Convert to Factor
researcher.data$blacks.less.than.deserve[researcher.data$blacks.less.than.deserve == -99] <- NA
researcher.data$blacks.less.than.deserve <- factor(researcher.data$blacks.less.than.deserve, 
                                                   labels=c("Strongly Agree", "Agree", "Disagree", "Strongly Disagree", "Neither Agree nor Disagree"))
researcher.data$blacks.less.than.deserve <- factor(researcher.data$blacks.less.than.deserve,levels(researcher.data$blacks.less.than.deserve)[c(1,2,5,3,4)])

researcher.data$blacks.like.other.minorities[researcher.data$blacks.like.other.minorities == -99] <- NA
researcher.data$blacks.like.other.minorities <- factor(researcher.data$blacks.like.other.minorities, 
                                                       labels=c("Strongly Agree", "Agree", "Disagree", "Strongly Disagree", "Neither Agree nor Disagree"))
researcher.data$blacks.like.other.minorities <- factor(researcher.data$blacks.like.other.minorities,levels(researcher.data$blacks.like.other.minorities)[c(4,3,5,2,1)])

researcher.data$blacks.not.trying[researcher.data$blacks.not.trying == -99] <- NA
researcher.data$blacks.not.trying <- factor(researcher.data$blacks.not.trying, 
                                            labels=c("Strongly Agree", "Agree", "Disagree", "Strongly Disagree", "Neither Agree nor Disagree"))
researcher.data$blacks.not.trying <- factor(researcher.data$blacks.not.trying,levels(researcher.data$blacks.not.trying)[c(4,3,5,2,1)])

researcher.data$blacks.slavery.discrimination[researcher.data$blacks.slavery.discrimination == -99] <- NA
researcher.data$blacks.slavery.discrimination <- factor(researcher.data$blacks.slavery.discrimination, 
                                            labels=c("Strongly Agree", "Agree", "Disagree", "Strongly Disagree", "Neither Agree nor Disagree"))

researcher.data$blacks.slavery.discrimination <- factor(researcher.data$blacks.slavery.discrimination,levels(researcher.data$blacks.slavery.discrimination)[c(1,2,5,3,4)])


researcher.data$racial.resentment <- ((as.numeric(researcher.data$blacks.less.than.deserve)-1)/4 + (as.numeric(researcher.data$blacks.like.other.minorities)-1)/4 +
                                      (as.numeric(researcher.data$blacks.not.trying)-1)/4 +  (as.numeric(researcher.data$blacks.slavery.discrimination)-1)/4)/4

### Model frame for predictions:
prediction.frame <- data.frame(treatment.black = c(1,0))

###### Estimate effect on racial resentment
resentment <- lm(racial.resentment  ~ treatment.black, data=researcher.data)
resentment_vcov <- vcovHC(resentment)
coeftest(resentment, resentment_vcov)
resentment_pval <- round(coeftest(resentment, resentment_vcov)[2,4], 3)
resentment_results <- data.frame(treatments = c(paste("Black Name\nn=",sum(model.frame(resentment)$treatment.black),sep=""), 
                                                paste("White Name\nn=",sum(model.frame(resentment)$treatment.black == 0),sep="")),
                                 point = predict(resentment, newdata=prediction.frame),
                                 se = sqrt(c(resentment_vcov[1,1] + resentment_vcov[2,2] + 2*resentment_vcov[1,2],resentment_vcov[1,1])))

resentment_results$lower95 <- resentment_results$point - qnorm(.975)*resentment_results$se
resentment_results$upper95 <- resentment_results$point + qnorm(.975)*resentment_results$se

#### Plot racial resentment
pdf("figures/FigureA8_Panel1.pdf", width=5, height=7.1)
p = ggplot(resentment_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="Racial Resentment Scale", limits=c(.2, .7))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",resentment_pval, sep=""))
p = p + ggtitle("Effect of researcher race\non racial resentment")
p = p + theme_bw1()
p
dev.off()

######## Black president
researcher.data$president.black.yes <- ifelse(researcher.data$president.black == 1, 1, 0)
researcher.data$president.black.no <- ifelse(researcher.data$president.black == 2, 1, 0)

black_president <- lm(president.black.yes ~ treatment.black, data=researcher.data)
black_president_vcov <- vcovHC(black_president)
coeftest(black_president, black_president_vcov)
black_president_pval <- round(coeftest(black_president, black_president_vcov)[2,4], 3)
prediction.frame <- data.frame(treatment.black = c(1,0))

black_president_results <- data.frame(treatments = c(paste("Black Name\nn=",sum(model.frame(black_president)$treatment.black),sep=""), 
                                                     paste("White Name\nn=",sum(model.frame(black_president)$treatment.black == 0),sep="")),
                                 point = predict(black_president, newdata=prediction.frame),
                                 se = sqrt(c(black_president_vcov[1,1] + black_president_vcov[2,2] + 2*black_president_vcov[1,2],black_president_vcov[1,1])))

black_president_results$lower95 <- black_president_results$point - qnorm(.975)*black_president_results$se
black_president_results$upper95 <- black_president_results$point + qnorm(.975)*black_president_results$se

#### Plot racial black_president
pdf("figures/FigureA8_Panel2.pdf", width=5, height=7.1)
p = ggplot(black_president_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% that would vote for a Black President", limits=c(.7, 1))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",black_president_pval, sep=""))
p = p + ggtitle("Effect of researcher race on\nwillingness to vote for\na black president")
p = p + theme_bw1()
p
dev.off()

####################

###### Women's Role
researcher.data$womenrole[researcher.data$womenrole == -99] <- NA
researcher.data$womenrole.equal <- ifelse(researcher.data$womenrole == 1, 1, 0)
researcher.data$womenrole.home <- ifelse(researcher.data$womenrole == 2, 1, 0)
researcher.data$womenrole.dk <- ifelse(researcher.data$womenrole == 3, 1, 0)

women_equal <- lm(womenrole.equal ~ treatment.woman, data=researcher.data)
women_equal_vcov <- vcovHC(women_equal)
coeftest(women_equal, women_equal_vcov)

women_equal_pval <- round(coeftest(women_equal, women_equal_vcov)[2,4], 3)
prediction.frame.woman <- data.frame(treatment.woman = c(1,0))

women_equal_results <- data.frame(treatments = c(paste("Female Name\nn=",sum(model.frame(women_equal)$treatment.woman),sep=""), 
                                                 paste("Male Name\nn=",sum(model.frame(women_equal)$treatment.woman == 0),sep="")),
                                      point = predict(women_equal, newdata=prediction.frame.woman),
                                      se = sqrt(c(women_equal_vcov[1,1] + women_equal_vcov[2,2] + 2*women_equal_vcov[1,2],women_equal_vcov[1,1])))

women_equal_results$lower95 <- women_equal_results$point - qnorm(.975)*women_equal_results$se
women_equal_results$upper95 <- women_equal_results$point + qnorm(.975)*women_equal_results$se

#### Plot women's equality
pdf("figures/FigureA9_Panel1.pdf", width=5, height=7.1)
p = ggplot(women_equal_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% that say women should be equal to men", limits=c(.7, 1))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",women_equal_pval, sep=""))
p = p + ggtitle("Effect of researcher gender on\nattitudes on gender equality\n")
p = p + theme_bw1()
p
dev.off()

######## Woman president
researcher.data$president.woman.yes <- ifelse(researcher.data$president.woman == 1, 1, 0)
researcher.data$president.woman.no <- ifelse(researcher.data$president.woman == 2, 1, 0)

woman_president <- lm(president.woman.yes ~ treatment.woman, data=researcher.data)
woman_president_vcov <- vcovHC(woman_president)
coeftest(woman_president, woman_president_vcov)


woman_president_pval <- round(coeftest(woman_president, woman_president_vcov)[2,4], 3)
prediction.frame.woman <- data.frame(treatment.woman = c(1,0))

woman_president_results <- data.frame(treatments = c(paste("Female Name\nn=",sum(model.frame(woman_president)$treatment.woman),sep=""), 
                                                     paste("Male Name\nn=",sum(model.frame(woman_president)$treatment.woman == 0),sep="")),
                                  point = predict(woman_president, newdata=prediction.frame.woman),
                                  se = sqrt(c(woman_president_vcov[1,1] + woman_president_vcov[2,2] + 2*woman_president_vcov[1,2],woman_president_vcov[1,1])))

woman_president_results$lower95 <- woman_president_results$point - qnorm(.975)*woman_president_results$se
woman_president_results$upper95 <- woman_president_results$point + qnorm(.975)*woman_president_results$se

#### Plot woman president - pilot
pdf("figures/FigureA9_Panel2.pdf", width=5, height=7.1)
p = ggplot(woman_president_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% that would vote for a woman president", limits=c(.7, 1))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",woman_president_pval, sep=""))
p = p + ggtitle("Effect of researcher gender on\nwillingness to vote for\na woman president")
p = p + theme_bw1()
p
dev.off()

#### Service spending
researcher.data$services.more <- ifelse(researcher.data$services == 2, 1, 0)
researcher.data$services.cut <- ifelse(researcher.data$services == 1, 1, 0)

service_cut_black <- lm(services.more ~ treatment.black, data=researcher.data)
service_cut_black_vcov <- vcovHC(service_cut_black)
coeftest(service_cut_black, service_cut_black_vcov)


service_cut_black_pval <- round(coeftest(service_cut_black, service_cut_black_vcov)[2,4], 3)
prediction.frame.black <- data.frame(treatment.black = c(1,0))

service_cut_black_results <- data.frame(treatments = c(paste("Black Name\nn=",sum(model.frame(service_cut_black)$treatment.black),sep=""), 
                                                     paste("White Name\nn=",sum(model.frame(service_cut_black)$treatment.black == 0),sep="")),
                                      point = predict(service_cut_black, newdata=prediction.frame.black),
                                      se = sqrt(c(service_cut_black_vcov[1,1] + service_cut_black_vcov[2,2] + 2*service_cut_black_vcov[1,2],service_cut_black_vcov[1,1])))

service_cut_black_results$lower95 <- service_cut_black_results$point - qnorm(.975)*service_cut_black_results$se
service_cut_black_results$upper95 <- service_cut_black_results$point + qnorm(.975)*service_cut_black_results$se

#### Plot service cut - black
pdf("figures/FigureA8_Panel3.pdf", width=5, height=7.1)
p = ggplot(service_cut_black_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% supporting greater service spending", limits=c(.3, .9))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",service_cut_black_pval, sep=""))
p = p + ggtitle("Effect of researcher race on\nwillingness to support\npublic service expansion")
p = p + theme_bw1()
p
dev.off()

### Gender effect

service_cut_woman <- lm(services.more ~ treatment.woman, data=researcher.data)
service_cut_woman_vcov <- vcovHC(service_cut_woman)
coeftest(service_cut_woman, service_cut_woman_vcov)


service_cut_woman_pval <- round(coeftest(service_cut_woman, service_cut_woman_vcov)[2,4], 3)
prediction.frame.woman <- data.frame(treatment.woman = c(1,0))

service_cut_woman_results <- data.frame(treatments = c(paste("Female Name\nn=",sum(model.frame(service_cut_woman)$treatment.woman),sep=""), 
                                                       paste("Male Name\nn=",sum(model.frame(service_cut_woman)$treatment.woman == 0),sep="")),
                                        point = predict(service_cut_woman, newdata=prediction.frame.woman),
                                        se = sqrt(c(service_cut_woman_vcov[1,1] + service_cut_woman_vcov[2,2] + 2*service_cut_woman_vcov[1,2],service_cut_woman_vcov[1,1])))

service_cut_woman_results$lower95 <- service_cut_woman_results$point - qnorm(.975)*service_cut_woman_results$se
service_cut_woman_results$upper95 <- service_cut_woman_results$point + qnorm(.975)*service_cut_woman_results$se

#### Plot service cut - woman
pdf("figures/FigureA9_Panel3.pdf", width=5, height=7.1)
p = ggplot(service_cut_woman_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% supporting greater service spending", limits=c(.3, .9))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",service_cut_woman_pval, sep=""))
p = p + ggtitle("Effect of researcher gender on\nwillingness to support\npublic service expansion")
p = p + theme_bw1()
p
dev.off()

#########################################################3
#### Does treatment affect finishing?

finish <- lm(NoFinish ~ treatment.whitewoman + treatment.blackwoman + treatment.blackman + treatment.whiteman, data=researcher.data.full)
finish_vcov <- vcovHC(finish)
coeftest(finish, finish_vcov)

prediction.frame.finish <- data.frame(treatment.whitewoman = c(1,0,0,0), treatment.blackwoman = c(0,1,0,0), treatment.blackman= c(0,0,1,0), treatment.whiteman = c(0,0,0,1))

finish_results <- data.frame(treatments = c(paste("White Woman\nn=",sum(model.frame(finish)$treatment.whitewoman==1),sep=""), 
                                            paste("Black Woman\nn=",sum(model.frame(finish)$treatment.blackwoman == 1),sep=""),
                                            paste("Black Man\nn=",sum(model.frame(finish)$treatment.blackman == 1),sep=""),
                                            paste("White Man\nn=",sum(model.frame(finish)$treatment.whiteman == 1),sep="")),
                             point = predict(finish, newdata=prediction.frame.finish),
                             se = sqrt(c(finish_vcov[1,1] + finish_vcov[2,2] + 2*finish_vcov[1,2],
                                         finish_vcov[1,1] + finish_vcov[3,3] + 2*finish_vcov[1,3],
                                         finish_vcov[1,1] + finish_vcov[4,4] + 2*finish_vcov[1,4],
                                         finish_vcov[1,1])))


finish_results$lower95 <- finish_results$point - qnorm(.975)*finish_results$se
finish_results$upper95 <- finish_results$point + qnorm(.975)*finish_results$se

#### Plot effect on completion
pdf("figures/FigureA10.pdf", width=7, height=7)
p = ggplot(finish_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% Failing to Complete Survey", limits=c(0, .2))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions", sep=""))
p = p + ggtitle("Effect of treatments on survey non-completion rates")
p = p + theme_bw1()
p
dev.off()

#####################################
############ Attention Checks

ac_pass <- lm(ac.pass.both ~ treatment.whiteman, data=researcher.data)
ac_pass_vcov <- vcovHC(ac_pass)
coeftest(ac_pass, ac_pass_vcov)


ac_pass_pval <- round(coeftest(ac_pass, ac_pass_vcov)[2,4], 3)
prediction.frame.ac <- data.frame(treatment.whiteman = c(1,0))

ac_pass_results <- data.frame(treatments = c(paste("White, Male Name\nn=",sum(model.frame(ac_pass)$treatment.whiteman),sep=""), 
                                             paste("Other Name\nn=",sum(model.frame(ac_pass)$treatment.whiteman == 0),sep="")),
                                      point = predict(ac_pass, newdata=prediction.frame.ac),
                                      se = sqrt(c(ac_pass_vcov[1,1] + ac_pass_vcov[2,2] + 2*ac_pass_vcov[1,2], ac_pass_vcov[1,1])))

ac_pass_results$lower95 <- ac_pass_results$point - qnorm(.975)*ac_pass_results$se
ac_pass_results$upper95 <- ac_pass_results$point + qnorm(.975)*ac_pass_results$se

#### Plot attention checks
pdf("figures/FigureA11.pdf", width=5, height=7)
p = ggplot(ac_pass_results, aes(y=point, x=treatments))
p = p + scale_y_continuous(name="% that successfully completed both attention checks", limits=c(.7, 1))
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_x_discrete(name=paste("Treatment conditions\np-value of difference = ",ac_pass_pval, sep=""))
p = p + ggtitle("Effect of researcher gender\n and race on attention\n check completion")
p = p + theme_bw1()
p
dev.off()
